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We present a study of the phononic thermal conductivity of isotopically disordered carbon nan- 
otubes. In particular, the behavior of the thermal conductivity as a function of the system length is 
investigated, using Green's function techniques to compute the transmission across the system. The 
method is implemented using linear scaling algorithms, which allows us to reach systems of lengths 
up to L = 2.5 fim. (with up to 400,000 atoms). As for ID systems, it is observed that the conductiv- 
ity diverges with the system size L. We also observe a dramatic decrease of the thermal conductance 
for systems of experimental sizes (roughly 80% at room temperature for L — 2.5 ^m), when a large 
fraction of isotopic disorder is introduced. The results obtained with Green's function techniques 
are compared to results obtained with a Boltzmann description of thermal transport. There is a 
good agreement between both approaches for systems of experimental sizes, even in presence of 
Anderson localization. This is particularly interesting since the computation of the transmission 
using Boltzmann's equation is much less computationally expensive, so that larger systems may be 
studied with this method. 

PACS numbers: 65.80.-|-n, 61.46. Fg, 63.22.-m, 63.20.kp 



I. INTRODUCTION 

Carbon nanotubes (CNTs) are very interesting mate- 
rials for nanoscale electronic devices due to their out- 
standing mechanical and electrical (depending on their 
chirality, CNTs can be semiconducting or metallic) prop- 
erties. Recently, it was also discovered that CNTs have 
very good thermal properties, as measured experimen- 
tally for individual single walled carbon nanotubes^ or 
as estimated by computer simulations (see the references 
in Section [il Bp . At room temperature, the thermal con- 
ductance of carbon nanotubes seems to be dominated 
by the phononic contribution, even for metallic carbon 
nanotubes.— 

Thermal properties are usually investigated using 
Fourier's law. For nonequilibrium steady states where 
the system is put in contact with two reservoirs at dif- 
ferent temperatures, there is a net energy flow from the 
hotter to the colder reservoir. The heat current density 
J is proportional to the temperature gradient 
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K being the thermal conductivity (a tensor, in general). 
Denoting by AT the temperature difference between the 
reservoirs, and by L the system size in the direction of 
the temperature gradient, 
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provided the temperature profile is linear. For usual three 
dimensional materials, the thermal conductivity does not 
depend on the system size, and so, it is a well-defined 
thcrmodynamical quantity. The situation is different for 



one-dimensional (ID) systems, for which the thermal flux 
is related to the transmission function of phonons from 
one reservoir to the other. For defect free periodic one- 
dimensional harmonic systems, there is no phonon scat- 
tering mechanism. Those systems can sustain a current 
which does not depend on the system's length,— so that 
|J|/AT is constant. The thermal conductivity therefore 
diverges as L, and is not well-defined. In general, one 
dimensional (ID) systems in which scattering processes 
can take place should exhibit an intermediate scaling 
I J I/I VT I ^ with < a < 1, in which case the thermal 
conductivity is again not well-defined. 

It is believed that CNTs should have a behavior rem- 
iniscent of ID systems, although such claims should be 
backed up by more systematic studies. The dependence 
of the CNT conductivities on the system length should 
therefore be a major and primary concern in any study 
of the thermal conductance. Only very few experimental 
studies on the length dependence of the thermal conduc- 
tance of carbon nanotubes are available to our knowl- 
edge (see for instance Rcf. |^ and Ref. 0), and numerical 
results are still rare (see the references in Section IIIBI) . 
More importantly, Fourier's law, which is not valid in 
general for ID systems, is sometimes assumed to hold to 
interpret experimental measurements in order to extract 
conductivities ^ 

In order to have a well-defined conductivity, a nec- 
essary condition (which may not be sufficient) is that 
some scattering mechanisms can take place, so that the 
transmission, hence the thermal flux, is reduced. Several 
scattering mechanisms exist in actual materials, which 
may be intrinsic, as for inelastic anharmonic phonon- 
phonon scattering; or extrinsic, as for elastic phonon scat- 
tering with defects. Experimental results showed that 
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CNTs of lengths 2.76 /xm may exhibit nearly ballistic 
transport ii This justifies that anharmonic scattering may 
be neglected if the elastic scattering processes with de- 
fects are important. The most simple defect that can 
be experimentally controlled is isotope disorder, which 
amounts here to replacing the usual ^^C atom by one its 
"'^^C isotope. CNTs with isotope disorder have already 
been synthetized^ and experimental results on boron ni- 
tride tubes^ showed that isotope disorder could lead to 
dramatic changes in the thermal conductivity. More- 
over, isotopically disordered harmonic systems are also 
the simplest systems that can be treated exactly with 
quantum statistics. 

There are several methods to compute the thermal con- 
ductance of isotopically disordered harmonic systems us- 
ing quantum statistics (we therefore disregard molecular 
dynamics techniques which give only results within the 
classical framework). These systems exhibit Anderson 
localization^ which arises from an interference effect 
of different scattering paths, and is thus a manifesta- 
tion of the ondulatory nature of the phonon vibrations. 
The Green's function technique is then a very appealing 
method to compute the thermal conductance in those 
systems, since it gives the exact transmission, has a com- 
putational cost scaling linearly with the system length, 
and is straightforward to parallelize since transport is 
coherent. This allows to reach systems of lengths up to 
L = 2.5 /im (with up to 400,000 atoms). A Boltzmann 
approach, which describes the evolution of the phonon 
distribution and therefore treats phonons as particles and 
not as waves, gives only an average description of the 
phonon flows in the system. It cannot account for the 
Anderson localization of states, and it is therefore un- 
clear whether it can provide a good approximation to 
the exact transmission. A Boltzmann treatment of the 
transmission is however very interesting since the method 
is computationally less expensive (no averages over the 
disorder realizations have to be taken, the computational 
scaling with respect to the CNT index is more favorable, 
and the inclusion of anharmonic scattering processes is 
easier than for Green's functions). 

In this paper, we will address three problems, system- 
atically studying the length dependence of the thermal 
conductivity with Green's function techniques: 

(i) does a well-defined conductivity (independent of the 
length) exist for harmonic disordered CNTs, or is 
the thermal transport anomalous? The theoretical 
results on heat transport in ID chains suggest that 
the transport is anomalous, but can the asymptotic 
divergence exponent be estimated, or should incred- 
ibly long tubes be considered before such a regime 
is found? 

(ii) how much is the thermal conductance reduced by 
isotopic disorder? How does this decrease depend 
on the temperature? 

(iii) in order to reach even larger systems (larger diame- 
ters and/or lengths), a third questions is whether 



the Boltzmann approach to thermal transport a 
good approximation for CNTs of experimental 
lengths! 

The paper is organized as follows. We present briefiy 
the notations and the numerical method in Section IIIII 
(see Appendices |X] and [B] for more precisions) . Wc then 
turn to isotopically disordered one dimensional chains in 
Section ITVl in order to test the validity of the Boltzmann 
approach in this simple case. Section |V] presents our new 
results on the systematic study of the length dependence 
of the thermal conductivity for CNTs. A conclusion is 
presented in Section IVTl 

II. HEAT TRANSPORT IN (QUASI) 
ONE-DIMENSIONAL SYSTEMS: A BRIEF 
REVIEW 

A. General features of heat transport in 
one-dimensional systems 

The theoretical and numerical results of heat transport 
in ID chains are not always acknowledged in computa- 
tional studies of thermal conductivities of CNTs. There 
have been many studies on the (non) validity of Fourier's 
law in one dimensional chains. There are two impor- 
tant review papers on this topic, written either from a 
mathematicalii or a physicali^ viewpoint. The main 
findings to this date may be summarized as follows: 

(i) there are chains in which Fourier's law can be shown 
to hold, in the sense that the conductivity does not 
depend on the system size. These situations arise 
when translational symmetry is broken (through the 
addition of some on-site potential in an anharmonic 
system for instance, or when momentum is not con- 
served by the dynamics, see the discussion in Ref.fisl 
for more precisions), or for specific models such as 
systems with self-consistent stochastic reservoirs at 
each lattice sites^ (a model recently studied math- 
ematically in Ref . [l5l ) ; 

(ii) there are chains with anomalous conductivities, i.e 
depending on the system size L. If the interaction 
between the particles on the chain depends only on 
their relative distances (and no on-site potential is 
considered, in particular), then anomalous heat con- 
duction is generally observed ( J/VT ~ L" diverges 
as some power-law), even with anharmonic inter- 
actions and/or mass disorder (see the references in 
Ref. [T^ . The precise exponent for systems with 
anharmonic interactions is still a matter of debate, 
with theoretical result o^^'^^ ranging from a = 1/3 
to a = 2/5, the most recent numerical results for 
Fermi-Pasta- Ulam chains giving values closer to 1/3 
(see Ref. Hi); 

(iii) in some models, the scaling of the conductivity de- 
pends on the details of the thermal reservoirs, the 
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way they are coupled to the system, or the bound- 
ary conditions usedj ^^i^^ 

Most of these results are found from numerical (molec- 
ular dynamics) simulations, though sometimes it is pos- 
sible to prove mathematically the (non) existence of an 
i-independent conductivity for the model at hand. 

More precise results on heat transport in isotopically 
disordered harmonic lattices (which is the focus here) 
are recalled in Section IIV Al when presenting numerical 
results for ID chains. 



B. Heat transport in CNTs 

The most popular method to compute conductivities is 
nowadays molecular dynamics (MD). A noticeable excep- 
tion is the Green's function approach to defect scattering 
of Refs. and S Sometimes, the Boltzmann equation 
may be used as well. We now give a brief review of MD 
studies. 

Within the MD framework, the most popular ap- 
proaches for thermal transport are (i) equilibrium meth- 
ods, which rely on the Green-Kubo formula; and (ii) 
nonequilibrium steady state (NESS) methods, where a 
non-zero heat flux is created in the system, either through 
boundary forcing (temperature difference) or using some 
mechanical forcingf^i Although the numerical procedures 
are by now standard computational methods in molecular 
dynamics, they are not free from the general limitations 
and issues concerning the convergence of the simulation 
results. For instance, it is known that the numerical com- 
putation of the conductivity relying on the Green-Kubo 
formula requires very long time integrals of the heat cur- 
rent autocorrelation function, and averages over many 
trajectories starting from different initial conditions. Be- 
sides, as the system size increases, the typical correlation 
time increases as welli^ In general, we expect a compu- 
tational cost scaling as O(L^) since one force iteration 
of the MD solver requires at least 0(L) operations, and 
the typical time (be it a correlation time in the Green- 
Kubo case, or the relaxation time to reach a steady state 
for NESS computations) is expected to scale as 0{L) as 
well. Moreover, to obtain the temperature dependence 
of the thermal conductivity, many simulations should be 
performed, each one at a different temperature. 

Critical reviews of MD results have been writteni^i^ 
In particular, the upper limit to the conductance given 
by the thermal ballistic conductance can be computed, 
and results of MD simulations compared to this upper 
limit »22, It is important to notice that MD results cannot 
be above the ballistic upper limit to the conductance. 
When this is the case, it is likely that MD results are not 
converged. Besides, Ref. [13 showed that in some cases, 
the simulation results depend on the boundary conditions 
and on the simulation technique used. 

In view of the MD computational scalings, only short 
nanotubes are computed (currently, the computations 
are usually restricted to system sizes L 0.2 /im. 



whereas experimental lengths are of the order L ~ 2 fim 
- however, some MD results for 3.2 /j,m CNTs are 
reported^^) and only few works studied the size depen- 
dence of the thermal conductivity for defect free CNTs 
in the framework of classical MD,^i^i25^i27i28 This 
amounts to considering only the effect of anharmonic- 
ity on the thermal conductance. The results, in agree- 
ment with the predictions of the one-dimensional heat 
transport theory, confirm the divergence of the conduc- 
tivity with increasing system sizes. More recently, the 
effect of defects and impurities on the thermal conduc- 
tivity have been investigated. This includes isotopic 
disorderi2^i22i22ii2£ and functionalizatioui^ In those cases, 
no systematic study of length dependence was performed. 



III. MODELS FOR THERMAL TRANSPORT 

We present in this section two important models of 
transport: (i) the Green's function approach, which 
solves the model in its full atomistic details, and is ex- 
act for harmonic systems; (ii) the computationally less 
expensive but more approximate Boltzmann treatment. 
Both approaches will be used here to compute the de- 
crease of the thermal conductance arising from a disor- 
dered region embedded in an otherwise perfect medium. 
We remark that the Green's function treatment gives the 
exact transmission for a given realization of the mass dis- 
order. Therefore, averages over the disorder realizations 
have to be performed, unless the transmission function 
does not vary much from one disorder realization to the 
other. On the contrary, the Boltzmann approach already 
gives some average transmission. 



A. The Green's function treatment of thermal 
transport 

1. Description of layered one-dimensional systems 

Coherent thermal transport, presented in a pedagog- 
ical fashion in Ref. [s^ . involves expressions very simi- 
lar to the usual expressions for electronic transport j^i^i 
The reason why Green's functions are introduced can 
be understood by explicitly integrating the equations of 
motions of the Heisenberg operators, projecting out the 
reservoir degrees of freedom, and turning to Fourier vari- 
ables since only the steady state of the system is of in- 
terest. Such a derivation follows the lines of the quan- 
tum Langevin equation,— and the relationship with the 
nonequilibrium Green's function (NEGF) formalism has 
been precised in Ref. [s^. Rigorous proofs of the existence 
of the NESS can be performed within the C* scattering 
formalismj2I 

The lattice thermal conductance arises from 
phononSf^i^ which are quantized displacements of 
harmonic lattices. We consider an infinitely extended 
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FIG. 1: View of a layered ID system. The infinite system 
is decomposed into three regions: a semi-infinite left lead, a 
semi-infinite right lead (which are both assumed to be per- 
fect), and a finite central disordered region. Assuming here 
that the atoms in the periodic unit cell interact with atoms 
located in the three neighboring cells on each side, a funda- 
mental "supercell", called a layer, can be considered. Atoms 
in such a layer interact only with atoms in the two adjacent 
layers. 



system as depicted in Fig. [T] The Hamiltonian is 

^(q,p) = ^q*^q + ^p*M-ip, 

where the infinite dimensional vectors 



q= (•■•,q»,i 



1 qi.wi qi+i,i5 



> Pi,JV, Pi+l,li 



stand respectively for displacements and momenta. The 
first index refers to the cell to which the atom belongs, 
and the second indexes the atom within a cell. The cell 
in question can be the periodic unit cell used to generate 
the system, but it may also be some convenient supercell 
(see below). The infinite dimensional matrix M is the 
(diagonal) mass matrix of the system: 



M 







The matrix K is the interatomic force constant matrix. 
It is assumed to be short ranged in the sequel, so that an 
atom interacts only with atoms located in a few neigh- 
boring unit cells. Changing coordinates to mass weighted 
coordinates, the transport properties of the system can in 
fact be completely characterized by the harmonic matrix 



(2) 



We restrict ourselves in this study to a disordered re- 
gion embedded in a perfect medium (see Fig. [1]). The 
case of mass disorder is then dealt with by considering the 
mass of the particles to be randomly distributed. In the 



most physical case, namely isotopic disorder, the proba- 
bility to have the mass m at a given site is 1 — c, and the 
probability to have a mass m -t- Am is c, where < c < 1 
denotes the disorder concentration. The assumptions on 
the system imply that the mass matrix is of the general 
form 



Mr 







A/,y, 



(3) 







M, 



R 



where in AI^ and Alfi all the masses are equal to m, while 
in Mgys they are randomly distributed. 

Carbon nanotubes are quasi one dimensional systems, 
that is, systems of finite size in the transverse directions 
and (infinitely) extended in the remaining direction. We 
consider a fundamental supercell (a layer) composed of 
possibly several unit cells. The number of cells in this 
fundamental structure is determined by the condition 
that an atom in a layer interacts only with atoms in the 
two adjacent layers. The harmonic matrix therefore has 
the generic block tridiagonal shape 



A = 



Al 1l 







(4) 







^fl Ar 



The (infinite) subsystems associated with the semi- 
infinite matrices A^, A^ represent some reservoirs to 
which the (finite) system Asys is coupled through the in- 
teraction terms Tl, T^. The block tridiagonal structure 
can be read off the expressions of the matrices (for sim- 
plicity, we assume that there is no mass disorder in the 
first and the last layer of the central region): 



At. = 







r a 
r* 
... 
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-2 Tn,^ -I 
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In the above expressions, a and r are 3iVat x 3A^at real 
matrices, A^at denoting the number of atoms in a layer. 
The matrix a represents the interactions within a layer, 
and the matrices r* , t model the interactions with respec- 
tively the left and right neighboring layers (For a matrix 
M., M.^ denotes the transpose matrix, while A^t denotes 
the hcrmitian conjugate of Ai in the sequel). If Mayers 
is the number of layers composing the disordered region, 
there arc A^at-^iaycrs atoms in the central part, and the 
matrix A^y^ is of size SA^^at^^iaycrs x 3 A^at Mayers- The ma- 
trices with underscripts ai,Ti, r* differ from the matrices 
a, r, r* because of the mass disorder (see Eqs. ([l])-©)- 



2. Heat current in terms of Green's functions 

The Green's function of the whole system is defined as 
the limit 

G+{lu) = lim {oj"^ + iri - A)-^ 

when this limit exists 4^ In numerical computations, the 
parameter is a small positive number. Knowledge of 
the Green's function implies the knowledge of the eigen- 
values and density of states of the system. An interesting 
feature of the Green's function formalism in the harmonic 
case is that the reservoir degrees of freedom can actually 
be projected out thanks to the linearity of the interac- 
tions, so that some effective dynamics only in the region 
of interest is recovered. The effective Green's function 
for the central region is 

G+M = lim + i'? - ^y. - ^tH - ^U^))"' , 

(5) 

where 

^+{uj) = lim T^(cc;2+i7y- 

j:+{iu) = lim ZR{uj^+ir]-Ai,)-H%. (6) 

The operators EJ, EJ are self-energies, which model the 
coupling of the system with the contacts. In particular, 
the imaginary part of the self energy is usually associated 
with some lifetime (due to phonons flowing out of the 
central region) . This is maybe more easily understood in 
the quantum Langevin framework, where the self energy 
is the Fourier transform of the friction kerneli^ 

When there arc no incoherent scattering processes, the 
current is given as the superposition of phonons going 
from the left reservoir to the right one, minus the flow 
of phonons going from the right reservoir to the left one. 
Therefore, phonons entering the system can be traced 
back to the reservoir where they come from and so, the 
summation is restricted to phonons with positive group 
velocity for phonons coming out of the left reservoir, and 
phonons with negative velocities for the ones coming out 
of the right reservoir. When the reservoirs are at different 



temperatures (respectively, and Tr), the heat current 
flowing through the system i o32,33,36 

ATl, Tn) = / ^r(^)(M H - hA^)) d^. (7) 

where the transmission factor T{lj) is 

r{Lo) = Tr [r+(a;)G+,(a;)r+(c.) (G+,(c.))t] . (8) 

In the above expressions, 

r+(^) = -2Im(E+(^)) {a^L,R) 

is related to the imaginary part of the self-energies and 
the functions /t are the Bose-Einstein distributions 



exp 



- 1 



The expression of the thermal conductance is exact since 
the system is harmonic. The practical computation of 
the transmission is presented in Appendix VK\ 

Introducing t(w) = ^L(w)l/2G'+,(^J)^J^(w)^/^ the 
transmission can be rewritten as T{ijj) = Tr(t([j)<(tj)t). 
It is therefore nonnegative. In fact, when there is no 
disorder, the transmission is ballistic and equal to the 
number of conducting channels at this pulsation. In the 
presence of mass disorder, the transmission is lower than 
the ballistic transmission. 



3. Thermal conductance and thermal conductivity 

The thermal conductance associated with the heat cur- 
rent ^ is obtained in the limit AT = Tr - ^ as 



9{T) 



lim 

AT^O 



J T 



AT 



.T 



AT 



AT 



27r' 



dT 



(9) 



A remarkable fact about the heat current is that the as- 
sociated thermal conductance is quantized;^ the quanta 
of thermal conductance being 



(10) 



This expression is obtained in the limit T — s- for a sin- 
gle acoustic branch in the ballistic regime. The quan- 
tum of thermal conductance has been experimentally 
measured4^ In the sequel, conductances will most often 
be normalized with respect to the quantum of thermal 
conductance: 



9{T) 



9{T) 

goiTY 



(11) 
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The dimcnsionlcss quantity g{T) is called "normalized 
conductance" in the results presented below. 

For ID systems, the thermal conductivity k is the ther- 
mal conductance divided by the system cross section A, 
and multiplied by the system length L: 



The notion of cross section depends on the system consid- 
ered and on the conventions used (see Section FVAI for our 
conventions in the case of CNTs). From this definition, 
it is however clear that the existence of a well-defined 
(length independent) non-vanishing thermal conductiv- 
ity requires a decrease of the thermal current (or thermal 
conductance) as 1/L. Since the thermal conductivity is 
not well-defined a priori, the thermal conductance is a 
much better notion to handle. 

The practical computation of the thermal conductance 
using ^ therefore boils down to computing the trans- 
mission factor T(tt') given by ([5|) for the whole phonon 
spectrum (see Appendix E]) . 

B. The Boltzmann approximation of coherent 
transport 

Boltzmann method is a less expensive but more ap- 
proximate method than the Green's function treatment 
for the computation of the thermal conductance. It al- 
lows to compute conductivities for very long systems, and 
therefore to study the length dependence of the conduc- 
tivity. For instance, the effect of anharmonicity (in terms 
of 3 and 4 phonon scattering processes) has been con- 
sidered for tubes of lengths up to L = 1 m4^ Recent 
derivations of the Boltzmann equation from microscopic 
dynamics have been presented in Refs. [U andl45l in the 
so-called kinetic limit (vanishingly small mass disorder or 
anharmoniciticS; long time limit, interatomic spacing go- 
ing to zero). The proofs arc sometimes only formal, the 
central object in the overall procedure being the Wigner 
distribution of phonons. The Boltzmann equation there- 
fore describes some average behavior of the phonons. 

The central object describing this average behavior is 
the phonon distribution nj(ijj,x,t). The index j labels 
the phonon branch, u) is the phonon pulsation, while x is 
the position along the system. The quantity nj{u;,x,t) 
therefore counts the number of phonons of pulsation uj 
in the j-th branch at a point x at time t. The evolution 
of the phonon distribution is governed by the Boltzmann 
equation, which consists of two terms: 



• a transport term, which models the free flow of 
phonons through the harmonic system of reference; 

• a collision term, which models the rate at which 
phonons from one branch are scattered into another 
branch due to the defects and/or impurities in the 
actual system. 

Since we do not consider anharmonic effects in this study 
but only mass defects, the energy is conserved in the col- 
lision processes and the only scattering mechanism is the 
scattering from a phonon in a given branch to a phonon 
in another branch at the same energy. 

1. The Boltzmann equation 

For the system introduced in Section lTlI A li a phonon 
at a wavevector k is an eigenvector of the dynamical ma- 
trix of the perfect system: 

D{k) ^a + TG-'^^'' +T*e"'''°, 

where /q is the lattice parameter (determined by the 
length of the periodic unit cell) and k belongs to the 
first Brillouin zone [— tt/Zq, tt/Zq]- More precisely, the j- 
th phonon is Uj = (e~''''"'''Uj)„gz where uj is such that 

D{k)uj ^ ujj{k)'^ Uj. 

Notice that since a and r are real matrices and ujj is a 
real number, U* is also a phonon since it is an eigen- 
vector of the dynamical matrix D{—k) for the pulsation 
ujj{—k) = LOj{k). For a given pulsation tU > 0, the asso- 
ciated phononic branches j at this pulsation are all the 
solutions of the equation 

L0]{k) = UJ^i-k) = Cj2, < J < 3iVat, (13) 

for some k in the Brillouin zone. We denote by = 
N{lu) the number of branches at this pulsation, and by 
nj^a{ijj,x,t) the density of phonons of the j-th branch 
(j = 1, • ■ • , N{uj) upon reordering). The index ct = ±1 
labels the solutions (k) = lD"^ depending on the sign of 
the phononic velocity: cr = +1 corresponds to k points 
such that 

^jikj.+i) = w, vj = ^(fcj- +i) > 0, 

whereas the velocity is negative at points = — fcj^+i, 
and is actually equal to —Vj. 
The Boltzmann equation reads 



{dt+av,d,)n,A^,x,t)= ^-'"•"^'.^''^"'^ nj,^,,{u,x,t) - ^^^■^^^'^^ n.^Aio, x,t). (14) 
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The left-hand side of the above equation is the transport 
part of the Boltzmann equation. Notice also that lo does 
not enter explicitly the above equation. It can therefore 
be considered as a continuous parameter indexing a set 
of independent equations. 

The matrix W = [W^(j,cr),(j',a')]i<jJ'<JV, <t,o-'=±i de- 
scribes the interactions between the branches, i.e. the 
collision term. It is such that for all (j, a) ^ (j', cr'), 



VK(,-<,),(,v,') > 0, 



W(. 



The above conditions have physical interpretations: (i) 
the first condition (non-negativity of the matrix ele- 
ments) means that the scattering term models the de- 
cay of phonons of one branch into phonons of another 
branch; (ii) the second condition ensures the conserva- 
tion of the total phonon number; (iii) the last condition 
is some symmetry condition with respect to the prop- 
agation direction of phonons. The precise form of the 
scattering matrix depends on the problem at hand; see 
below the expression (fT5|) in the case of isotopically dis- 
ordered lattices. 



2. The collision term 

The expression of the scattering matrix are derived 
from a perturbative approach (based on the Fermi 
Golden Rulc)!i&^>i^ 



Wi 



Var(m) _2 
'0 , ,o a; 



El 

1=1 



ril) 



(15) 

where (m) and Var(m) = (m^) — (m)^ are the average 
mass and the variance of the mass disorder respectively. 
The three dimensional vector Uj^^(l) is the part of Uj 
corresponding to the three degrees of freedom of the l-th 
atom in the unit cell for phonon displacements Uj com- 
puted for a perfect lattice. 

The scattering term (fT5)) agrees with the mathematical 
results obtained by the kinetic limit of the microscopic 
dynamics in the case of a simple one-dimensional chain4^ 
From a numerical viewpoint, the expression (jlSp is very 
interesting since it allows an analytic computation of the 
scattering rates once the phonon spectrum has been com- 
puted. Appendix IB] presents the details of the numerical 
solution of (fTi|) using the scattering term ([TS]). 



3. Transmission properties using the Boltzmann equation 

Thermal properties are computed for nonequilibrium 
steady states where a heat current flows through the sys- 
tem. This is done by setting appropriate boundary con- 
ditions, and waiting for the system to equilibrate. The 



boundary conditions for are suited to the transmis- 
sion of a single phonon from the left (hot) reservoir to 
the right (cold) one: 

nj^+{LL>,0,t) = 1, nj__(a;, i, t) = 0. 

The numbers nj^-\^{uj^ L^t) and nj^-{uj^Q,t) are re- 
spectively the proportion of transmitted and reflected 
phonons. These proportions are computed using the 
Boltzmann equation, and this defines the transmission 
factor. 

When a stationary regime is reached {dtUj^a = for 
all J, (t), the phonon distributions do not depend on time 
anymore, and we drop the variable t in our notations. 
It is also easily seen that the transmission coefficient is 
independent of the position x, so that 



T{uj,x) = N{uj) 



N{uj) 



N{uj) 



N{uj) 



JV(w) 

3 = 1 

N(uj) 



(16) 



This coefficient is the Boltzmann equivalent of the trans- 
mission function T{uj) computed using Green's function 
techniques. 



IV. ONE DIMENSIONAL CHAINS 

Perfect onc-dimensional harmonic lattices with nearest 
neighbor interactions are described by the Hamiltonian 



H 



y- 



1 



2m 



k{q 



This model fits in the general framework presented in 
Section [IIII when taking iVat = 1 atoms per unit cell and 
layers of Mayors = 1 Unit cell. The intralayer interactions 
arc then described by the 1x1 matrix a = [2k /m], while 
the interlayer interactions t = = [— fc/m]. 

In this section, we will work in reduced units of masses 
and lengths, and consider a one dimensional chain with 
unit lattice spacing and particles of masses 1. Isotopic 
disorder is modelled by replacing with probability < 
c < 1 the mass of a particle by 1 + (5 with 5 ~ Am/m. 



Heat transport in isotopically disordered 
harmonic lattices 



Heat transport in isotopically disordered harmonic lat- 
tices has been thoroughly studied i^^i^^i^'-'i^^i^^ In partic- 
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ular, the case a finite disordered chain connected to two 
semi infinite perfect chains has been considered^ and 
the corresponding theoretical results were rederived from 
a different perspective and extended in Refs. [s^ and [TqI. 
It is suggested in Ref. \^ that the thermal conductiv- 
ity diverges as \/L [L being the chain length). This 
has been shown rigorouslj^ for an analogous continuum 
wave model (open boundary conditions and disorder). 

Heuristic arguments can be employed to back up the 
^/L divergence of the thermal conductivity. Denoting ex- 
plicitly the length dependence of the transmission func- 
tion by Tl{uj), the exact transmission is rewritten as 

rL(c^) ~exp(-L//ioc(c^)), (17) 
since it can be shown that the limit 

lim - y In Tl{u})= \ 

exists for all w > 0. This result was proved in Ref. [l^ 
(but earlier obtained in the limit of low concentration of 
defects^) . The proof is based on a theorem by Fursten- 
berg on the product of random matricesi^ The behavior 
of /loc around lo^ = can also be precised] For iso- 
topic disorder, it can be shown that 



lim 



1 



Var(7Ti) c(l — c) 



-♦0 LO- 



{1 + C5I2Y 



(18) 

This shows that high a frequency implies a small ^loc and 
therefore the associated eigenmodes are strongly local- 
ized, so that only the low frequency modes contribute to 
transport; in fact, only a fraction 0(L^^/^) of them since 



Var(m) , 
exp I , ; ^ Lu^ 



when uj ^ Q. This explains therefore the L 
the thermal conductance. 



'1/2 



decay of 



B. Comparison of the Green's function and 
Boltzmann treatments 



The transmission is computed using the Green's func- 
tion approach at A'^^^ = 5000 points uniformly spaced in 
the range [0,Wmax], with rj/uj^^^^ = 10~^^ (the self ener- 
gies having an analytical expression, see for instance^). 
The averages are taken over A^jisordcr = 100 realizations 
of the isotopic disorder for chains of length L = 10^, 

while iVdisordcr = 1000 for L = 10^ and A^disordcr = lO'^ 

for L = 10^ — lO"'. The average transmission functions 
computed from ^ arc presented in Fig. [2] in the case 
6 = 1/12 and c = 0.5. The mass variation S is the mass 
variation corresponding to substituting ^-^C with ^'^C. 

The transmission predicted by the Boltzmann ap- 
proach can be computed analytically in this simple case 
since there is only a single branch and there are there- 
fore only two conducting states for a given pulsation (due 
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FIG. 2: (color online) Averaged transmission function Tl{uj) 
as a function of uj, for different system sizes. From top to 
down: increasing system sizes from L = 100 to L = lO*". The 
black curves are the transmissions computed with a Green's 
function approach, and the red curves are obtained from the 
Boltzmann formula (1191). 



to the symmetry uj(k) 
Appendix [C| that 



u!{—k)). It can be shown (see 



Tl{u^) = 1 



'Boltz 



with 



1 Var(m) 
'^^^Boitz(w) {my 



(co) 



1 - 



(19) 



Notice that this expression of the transmission agrees at 
second order in uj with (fT7)l in the limit iv ^ 0. 

The comparison between the transmissions computed 
with the Green's function and the Boltzmann approaches 
shows that the Boltzmann treatment is a reasonable ap- 
proximation for low frequency modes, but predicts a 
transmission which is too large for higher frequencies. 
The critical frequency at which the Boltzmann transmis- 
sion starts to depart significantly from the Green's func- 
tion transmission decreases with the system size. There- 
fore, we expect the conductances to agree in the low tem- 
perature regime for all system lengths, and the relative 
error to increase with the system length at larger tem- 
peratures where the higher frequency part of the trans- 
mission are taken into account. 

This is indeed confirmed by the scalings of the normal- 
ized thermal conductance in Fig. [3l for different values 
of a; = fajj/k-QT. As expected, the asymptotic scaling 
L^^l"^ for the thermal conductance is reached for systems 
long enough when the transmission is computed using the 
Green's function approach. This could be anticipated 
since the Green's function results are exact for harmonic 
systems. On the other hand, the conductances predicted 
by the Boltzmann approach for systems long enough, are 
larger than the conductances obtained with the Green's 
function approach. The Boltzmann treatment does not 
predict the right asymptotic scaling either, though the 
discrepancy is not too large. This is due to the fact that 
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FIG. 3: (color online) Conductances computed using the 
Green's function approach and a Boltzmann treatment for 
two values of a; = TvjJ /ksT. The asymptotic regime where 
g r-j L~" is attained only for x = 1, with a ~ 0.48 for the 
Green's function approach, and a ~ 0.44 in the Boltzmann 
case. 



with the increasing size of the system (though the range 
of pulsations where the variance is not small decreases). 
However, the thermal conductance, which is the inte- 
gral of the transmission function, does not vary much 
from one realization of the mass disorder to the other. 
Fig. m (Top) presents a plot of a the transmission func- 
tion for a given realization of the mass disorder, which is 
compared to the average transmission, while Fig. 3] (Bot- 
tom) plots the variance of the transmission as a function 
of the pulsation. Therefore, it is not true that, when 
the size of the disordered region increases, all the trans- 
mission functions look similar (no thermodynamic limit 
seems to be reached) . There remains some intrinsic vari- 
ability, which can be used to determine exponential lo- 
calization IcngthsjS This is also reminiscent of the fact 
that, in molecular dynamics studies of harmonic disor- 
dered chains, non smooth steady temperature profiles are 
obtained when only a single realization of the disorder is 
consideredi^ 
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FIG. 4: (color online) Top: Transmission computed using 
the Green's function technique for a one dimensional chain of 
length L = 10^: average transmission (black) and transmis- 
sion for a single realization of the mass disorder (red). Bot- 
tom: Variance of the transmission functions computed using 
the Green's function techniques as a function of the pulsation, 
for different system lengths. 



the Boltzmann transmission does not decay fast enough 
with the system size for larger values of oj. 

An interesting feature of the transmission function in 
the one dimensional case is that the maximal variance 
of T{uj) over the realizations of the transmission com- 
puted with different mass disorders does not decrease 



V. NUMERICAL RESULTS FOR CARBON 
NANOTUBES 

A. Description of the model 

We use two distinct models to simulate CNTs. With 
the first one, labelled as 'ab-initio' in the sequel, the 
tube is simulated using the interatomic force con- 
stants (IFCs) computed from density functional the- 
ory calculations using the PWCSF package of the 
QUANTUM-ESPRESSO distributioii^ on a (5,5) CNT 
of experimental (curved) geometry.— With the second 
one, labelled as 'empirical' in the sequel, the tube is sim- 
ulated with a flat geometry using the empirical force con- 
stants of graphene given in Ref. 160 . The empirical model 
is used to simulate (5,5) and (10,10) CNTs. In any cases, 
the models are used consistently within both the Green's 
function and Boltzmann approaches. We restrict our- 
selves to armchair nanotubes for simplicity. Armchair 
nanotubes are very important for applications since they 
are metallic (here we compute only the phononic thermal 
conductance). We expect the results presented below to 
be robust with respect to the system chirality, because 
the ballistic conductance does not change much with the 
chirality 1^ 

In order to have a block tridiagonal form for the in- 
teraction matrix within the ab-initio model, only inter- 
actions within a cut-off radius i?cut = 12 Bohrs are taken 
into account. In this case, it is necessary to consider 
a supercell composed of two elementary cells depicted in 
Fig.[5]as a unit building block. The projection procedure 
of Ref. 1^ was used to ensure that acoustic sum rules are 
satisfied for the truncated IFC matrix. 

When the empirical model of Ref. [13 is used, carbon 
nanotubes are modeled as graphene nanoribbons with 
periodic boundary conditions in the transverse direction, 
and interactions up the third neighbors are taken into 
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FIG. 5: Armchair nanotube, with interactions up the third 
neighbor. The unit ceU of the periodic structure is represented 
by the box. 

account. Due to these simplifying assumptions, the in- 
plane and out-of-plane motions can be computed sepa- 
rately since their contributions arc independent. In this 
case, the elementary cell depicted in Fig.[5]can be chosen 
as a unit cell. Although less precise (though, as will be 
shown below, many results are in very good agreeement 
with the results obtained with the ab-initio based com- 
putations), this model is less expensive and was therefore 
used to study CNTs of larger diameters. 

We compare in Fig.[6]phonon dispersions for (5,5) arm- 
chair tubes obtained with the empirical model and with 
ab-initio computed IFC. The agreement is fair enough, 
though the number of acoustic modes is not correct, and 
the lowest modes do not have a quadratic dispersion re- 
lation as is the case for real nanotubes. Despite that, the 
thermal transport results for both models should agree 
quantitatively as shown by a comparison of ballistic con- 
ductances computed with the empirical model used here 
and ab-initio reference results. 

Actually, a more intrinsic property is the ballistic ther- 
mal conductance divided by the cross sectional surface of 
the material g/A (this is the so-called ballistic conductiv- 
ity per unit length, compare with Eq. (fT2|) ). The cross 
sectional area A is obtained from a "fattened" carbon 
ring: 

A = 2TTRd ^ 3nd, (20) 

where n is the index of the nanotube, R is the radius 
of the carbon ring and d = 3.35 A (the graphite inter- 
layer separation) is some characteristic length defining 
the width of the annular domain enclosing the carbon 
ring. Fig. [7] presents a plot of g/A as a function of the 
temperature, for two CNTs of different indices, as well 
as a reference curve computed from ab-initio results, and 
experimental resultsji 

Several conclusions can be drawn from this picture: 




k [Ti/a] 

FIG. 6: Phonon spectrum of a (5,5) CNT. Top: Ab-imtio 
computed interatomic forces. Bottom: Empirical model de- 
scribed in Section W Al 

(i) there is a good agreement between the empirical 
model used and the reference results computed with 
ab-initio interatomic force constants on the whole 
range of temperatures; 

(ii) the agreement with the experimental results from 
Ref. Ill for temperatures around room temperature 
and below, suggests that the present models con- 
tain the essentials ingredients for describing ther- 
mal transport. In particular, other effects (presently 
not taken into account) such as anharmonic interac- 
tions may be neglected in this temperature regime 
(and for lower temperatures). The measurements 
from Ref. [l] have been transformed into conductiv- 
ities per unit length by assuming that the CNTs 
used for the experiments have a diameter of 1 nm. 
The authors of Ref. [l] indeed precise that the CNTs 
used have a diameter in the range 1-2 nm, occa- 
sionally 2-3 nm. If the diameter is indeed 1 nm, 
then the experimental results are close to the ballis- 
tic conductance curve, which means that the ther- 
mal transport is nearly ballisticji If the diameters 
are larger, there is a reduction of about 50% of the 
conductance for the CNTs of experimental lengths 
(about L = 2.76 ^m) with respect to the ballistic 
conductance. This reduction may be attributed to 
anharmonic effects. We remark that, even in this 
case, the effect of isotopic disorder should still be 
noticeable since the results on the conductivities in 
Section |VE] show that isotope disorder can lead to a 
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FIG. 7: (color online) Ballistic conductivity per unit length 
(expressed in A), with the cross sectional area pO|l for a (5,5) 
(red) and a (lOJO) (black) nanotube. Experimental points 
taken from Ref. [l| are also reported. 



reduction of about 80% of the conductance at room 
temperature for CNTs of experimental lengths. 

In the sequel, results for temperatures much larger than 
room temperature are also presented. We warn however 
the reader that these results are given for completeness 
(having computed the transmission, the thermal conduc- 
tance can straightforwardly be evaluated at any temper- 
ature). It is likely that anharmonic effects become more 
important as the temperature increases, invalidating the 
use of a harmonic description in this temperature range. 



Parameters used 

The parameters used in this study are given in Table H] 
when the empirical model is used. In the case of ab- 
initio computed IFC, no averaging was performed (since, 
as shown below in Section IV D| the conductance does 
not change much from one realization of the mass disor- 
der to another) , but the CNTs considered have the same 
lengths. Since the interatomic distance at equilibrium is 
Lq ~ 1.44 A, the nanotubes considered have lengths up 
to L = 2.49 /im, which is a length typical of the cur- 
rently available experimental single walled CNTs. The 
larger system {i.e. (10,10) CNT of length L = 2.49 n^i) 
consists of 400, 000 atoms. 

In all this study, c = 0.5, and the mass disorder cor- 
responds to replacing ^^C by ^^C (Am/m = 1/12). The 
qualitative features presented are robust with respect 
to a smaller disorder concentration, but longer tubes 
should then be studied. Also, the regularizing param- 
eter r?/cj2iax = 10"® • 

For ab-initio computed IFC, the transmission is com- 
puted for Nf^ = 200 points uniformly spaced in the range 
< w < Wmax = 1650 cm-\ For (5,5) CNTs de- 
scribed by the empirical model, = 1000 points uni- 
formly spaced in the range < uj < LOmax ~ 1650 cm~^. 
For (10,10) CNTs described by the empirical model, the 
spectral region < w < 330 cm~^ is discrctized using 



TABLE I: Number of averages over the mass disorder real- 
izations as a function of nanotube length for CNTs described 
by the empirical model. 



Length (nm) 


2490 


823 


249 


82.3 


24.9 


8.23 


-^layers 


10* 


3300 


1000 


330 


100 


33 


^ disorder 


1 


3 


10 


30 


100 


300 



iV^ = 200 points, while N^^ = 200 other points are used 
for the remaining part of the spectrum. Such a decom- 
position ensures that the low frequency part of the spec- 
trum, which is very important to have reliable estimates 
of the thermal conductivity, is treated accurately, while 
keeping the overall computational cost reasonable. The 
Boltzmann transmissions are computed in all case with 
= 1000. 

We studied the dependence on the results on the dis- 
order realization, in the case of CNTs described by the 
empirical model. For each tube length L, A'disordcr real- 
izations of the isotopic disorder are considered, and the 
transmission arc averaged over the different realizations. 
We considered a fixed computational cost for each tube 
length, so that XA'disordor is constant. 



B. Average transmissions 

The average transmission functions are presented in 
Fig. [51 The different transmission pictures obtained show 
that the higher frequency modes are damped out first, 
while the acoustic modes are less perturbed. In general, 
the transmission decreases with the pulsation. Longer 
system sizes are necessary to modify substantially the 
transmission of acoustic modes. As the carbon nanotube 
index (equivalently, its diameter) increases, the spikes in 
the transmission become less marked, and the behavior of 
the transmission as a function of the pulsation is almost 
monotonic. 

Although not shown here, the maximal variance over 
the transmission functions computed for different realiza- 
tions of the mass disorder is constant for different tube 
lengths, as is the case for one-dimensional chains. How- 
ever, the relative variation of the transmission (defined 
as the square root of the variance divided by the aver- 
age transmission) decreases with the CNT index. This 
is due to the fact that transmission fluctuations remains 
of order 1, while the ballistic transmission (the number 
of conducting channels) increases proportionally to the 
system index. 

The Boltzmann transmission is in fair agreement with 
the Green's function transmission, especially at low fre- 
quencies. At higher frequencies, the Boltzmann trans- 
mission is usually higher than the Green's function one. 
This is in analogy with the results for the dimensional 
chains. The largest discrepancies arise however for pul- 
sations in the range cj = 400 — 500 cm~^, corresponding 
to the out-of-planc phonon modes. The agreement seems 
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FIG. 8: (color online) Transmissions obtained with the 
Green's function formalism (black curves) and the Boltzmann 
approach (red curves). In each of the three figures, the line 
corresponding to the largest transmission is the reference bal- 
listic transmission. The other lines correspond to tubes of 
different lengths, namely L = 25 nm, L = 249 nm and 
L = 2.49 /im, the highest transmissions corresponding to the 
shortest tubes. 



however better when ab-initio IFCs are used. This may 
be due to the fact that in-plane and out-of-plane modes 
interact in this case (indeed, the Boltzmann treatment is 
expected to be more accurate when the number of inter- 
acting modes is larger). 



C. Determination of the transmission regime at 
some chosen frequencies 

Anderson localization^ is often observed in one- 
dimensional disordered systems. In this case, the trans- 



mission function decreases exponentially with the sys- 
tem size (in average). This is the case for isotopi- 
cally disordered harmonic one-dimensional chain (see 
Section flV Ap . for which however the localization length 
goes to infinity when the pulsation goes to 0. When many 
conducting channels are present in a disordered harmonic 
system, it is unclear whether the transmission is still lo- 
calized, or if it becomes closer to a diffusive behavior. 

We discuss here whether the transmission regime is 
diffusive or localized for a given pulsation lo. There are 
several ways to determine the nature of the transmission 
regime, see for instance the recent work^ where diffu- 
sive and localized lengths for different system sizes are 
computed. We compute here the diffusive and local- 
ized lengths by resorting to some functional ansatz of 
the transmission function which allows to interpolate be- 
tween the two regimes. The transmission considered is 
the average transmission since this is the quantity of in- 
terest. The (average) transmission function computed by 
the Green's function and the Boltzmann approaches may 
be approximated by the following ansatz: 



'(w) = rbamstic(w) 1 + 



L 



^diff(w) 



(21) 

where Zdiff(w) is some characteristic length associated 
with diffusive transport at the pulsation uj (inducing a 
1/L decay) while l\oc{^) measures the propensity of the 
system to localize its states (inducing an exponential de- 
cay of the transmission as in the Anderson model). 

To determine the lengths Zdiff('^), ^ioc('-^)! we compute 
A^reaiizations valucs of thc transmission at a fixed value 
of uj for several lengths L ~ ii, . . . , iAfiongthaJ ^'^d P^'^" 
form a least-square fit of the average transmission based 
on the ansatz (|2ip . Many realizations should be con- 
sidered in order to have reliable results, and so, we re- 
stricted ourselves to the simple empirical model; here 
A^rcaiizations = 300. Bcsidcs, it is actually not always 
sufficient to perform a least-square fit directly on the 
transmissions. Indeed, when the transmission regime is 
(close to) localized, the direct estimates give poor results 
since the transmission is asymptotically very small, and 
the corresponding points have very small weights in the 
direct least-square fitting procedure. Therefore, we re- 
sort to a nonfinear least-square fit, where /(Ti(a')) and 
j^yansatz^^-j-j comparcd. The function / puts more 
emphasis on the asymptotic regime, i.e. on the smaller 
values of the transmission. The choice / = — log is con- 
venient, but other choices such as f{x) = \/x yielded 
similar results. 

Fig. [5] presents a typical fit in thc diffusive regime at 
the frequency w = 716 cm~^ for a (10,10) CNT, and com- 
pares thc normalized transmissions Ti(w)/Tbaiiistic('^) 
computed using both Boltzmann and Green's function 
methods. As can be seen, the functional form of thc 
ansatz (PT|) is indeed appropriate. 

Thc characteristic lengths as for different pulsations 
arc presented in TableHHfor (5,5) CNTs, and in Table|TlT] 
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FIG. 9: (color online) Fit of — ln(rL/rbaiiistic) as a function 
oiLatLu ^ 716 cm"^ for a (10,10) CNT described by the em- 
pirical model. The transmission computed with the Green's 
function and Boltzmann methods can be fitted with a good 
precision using the ansatz (|2ip (solid black and red lines). 



for (10,10) CNTs (for some pulsations, the transmission 
decreases too fast with increasing system size, and no 
fitting could be performed). Using these results, it is 
possible to determine whether the transmission regime 
is diffusive or localized. Notice however that the diffu- 
sive and localized lengths as defined by (PT|) should not 
be compared directly when they are of the same order 
of magnitude. Rather, given some attenuation factor 
< 7 < 1 (typically, 7 — 0.1), one should compare 
the lengths of the sample required to reduce the ballistic 
transmission by a factor 7. For the diffusive term, this 
is L ~ (1/7 — l)Zdiff(w), while for the locaHzed regime 
L ^ -Zioc(t^) 1117. 

The results of Tables [II] and IIIII show that, for the 
exact (Green's function) transmission, the localization 
length decreases with increasing frequencies, so that lo- 
calization effects become more and more important as 
the pulsation increases. This explains the fast decrease 
of the high frequency modes. The behavior of the diffu- 
sive lengths is not as clear cut as the behavior of the 
localization lengths: the diffusive lengths decrease at 
smaller pulsations, and then settle down to values around 
100 — 300 nm, which are therefore completely relevant for 
CNTs of experimental sizes. It seems that some minimal 
diffusive length is obtained when the number of channel 
is maximal, which is consistent with the picture that dif- 
fusive transport arises from random scattering between 
many channels. 

For the Boltzmann transmission, only the diffusive 
part is relevant since the localization lengths computed 
are always much larger than the diffusive lengths. There 
is also a trend towards shorter diffusive lengths as the 
pulsation increases. This may be explained by the fact 
that the scattering for higher frequencies is much more 
efficient than for acoustic modes, which, in turn, comes 
from the dependence in the expression of the scatter- 
ing matrix (fT5|) . 

Finally, it is observed that both diffusive and localiza- 
tion lengths increase with increasing diameter of the sys- 



TABLE II: Localization and diffusive lengths (expressed in 
nm) for a (5,5) CNT described by the empirical model. 
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TABLE III: Localization and diffusive lengths (expressed in 
nm) for a (10,10) CNT described by the empirical model. 



Green Boltzmann 
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tem, especially the localization length. This can related 
again to an increase of the number of available modes 
with increasing diameters, and thus, a reduced scatter- 
ing. 

D. Asymptotic scaling of the thermal conductance 

Before discussing the properties of the thermal conduc- 
tance, we first check whether the realizations of the mass 
disorder have a noticeable influence on the value of the 
normalized thermal conductance g given by (fTTj) . This 
is important from an experimental viewpoint since only 
few (if not only one) realizations of the mass disorder can 
be considered. Fig. [10] presents the normalized standard 
deviation of the conductance, for different system sizes, 
as a function of the temperature. The variations of the 
conductance are lower than 1% from one realization to 
the other. This is due to the fact that the conductance is 
the integral of the transmission function (with the proper 
weighting function), and therefore the variations of the 
transmission function for only one realization of the mass 
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FIG. 10: (color online) Standard deviation of the conduc- 
tance a{g) of a (5,5) CNT (described by the empirical model) 
divided by average conductance (g) , as a function of temper- 
ature, for different system sizes. 



Fig. [TT] presents the averaged normalized thermal con- 
ductance as a function of the system length (in log-log 
scale) for (5,5) and (10,10) CNTs. Three temperatures 
are considered in both cases: T = 50 K, T = 300 K, 
T = 1000 K. The scahngs obtained for T = 50 K show 
that the transmission regime is quasi-ballistic since the 
conductance does not change much with the system size. 
Indeed, at low temperatures, only the low frequencies 
modes, which are almost unaffected by the disorder, mat- 
ter. For higher temperatures, more and more higher 
frequency modes are introduced, so that disorder has 
a noticeable effect, and the conductance decreases as a 
function of the system size. When this makes sense, the 
slope a of the different curves have been estimated us- 
ing some least-square fitting, in order to characterize a 
power law decay g ~ The exponents found are in 

the range a = 0.4 — 0.5 at T = 300 K. The associated 
thermal conductivities therefore have a power-law diver- 
gence gL/A ^ with (3 = 1 — a = 0.5 — 0.6 for the 
range of lengths considered in this study. The general 
trends in the curves are quite comparable for (5,5) and 
(10,10) CNTs, independently of the model used. 

Notice that the asymptotic regime is not yet attained 
for nanotubes of lengths up to 2.5 /im (which are typi- 
cal experimental lengths). In this regime the exponent 
a is not universal since it depends on the tube length 
and on the amount of mass disorder. As in the case of 
one dimensional chains, the truly asymptotic regime cor- 
responds io 'g ^ (or k ~ i^^^), with associated 
transmission profiles where only acoustic modes have a 
non-zero transmission. To this end, much larger system 
sizes (or a larger mass disorder) should be considered. 
In any cases, it is observed that there is no well-defined 
conductivity. 
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FIG. 11: (color online) Variation of the normalized conduc- 
tance for temperatures T = 50 K (black curves), T — 300 K 
(red curves), T = 1000 K (blue curves), for the conductivities 
computed with the Green's function approach and the con- 
ductivities computed with the Boltzmann approach. For a 
given temperature, the Boltzmann curve is always above the 
Green's function curve. Top: (5,5) CNT, estimated a — 0.40 
at T = 300 K and a = 0.50 at T = 1000 K. Middle: (10,10) 
CNT, estimated a = 0.41 at T = 300 K and q = 0.50 at 
T = 1000 K. Bottom: (5,5) CNT with ah-imtio IPCs, esti- 
mated a = 0.43 at T = 300 K and a = 0.55 at T = 1000 K. 



E. Thermal conductance as a function of the 
temperature 

The temperature dependence of the normalized con- 
ductance is presented in Fig. [T^l Those curves are ob- 
tained for the transmission computed with the Green's 
function approach, and the curves for the transmission 
computed using a Boltzmann treatment have a very sim- 
ilar behavior with respect to temperature. As expected. 
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the ratio of the thermal conductance to the baUistic one 
is always smaller than 1, and converges to 1 in the low 
temperature limit. 

Several conclusions can be drawn from this picture: (i) 
the relative reduction of the conductance due to the iso- 
tope disorder does not depend on the CNT diameter; (ii) 
the decrease predicted by Boltzmann equation is in excel- 
lent agreement with the decrease found with the Green's 
function method; (iii) isotopic disorder can be very effi- 
cient in reducing the thermal conductivity, especially for 
tubes of experimental lengths, even at moderately high 
temperatures. For instance, for (5,5) and (10,10) tubes of 
experimental lengths {L = 2.49 /xm), the thermal conduc- 
tance is decreased by 80% at room temperature. Notice 
that the decrease in the thermal conductivity increases 
with the temperature. This is consistent with the re- 
sults of the previous sections since, as the temperature 
is increased, more and more higher frequency modes are 
introduced, and the transmission function is almost a de- 
creasing function of the pulsation. 
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FIG. 12: (color online) Ratio of the thermal conductance to 
the ballistic conductance, for (10,10) empirical CNT (blue 
curves), (5,5) empirical CNT (red curves) and (5,5) ab-initio 
CNT (black curves). From top to down: increasing tube 
lengths. 

The thermal conductivity (or conductance) predicted 
by the Boltzmann approach is compared to the conduc- 
tivity (or conductance) obtained by the Green's function 
approach. The results are presented in Fig. 1131 where the 
error 



e(T) 



KBoltz 



is plotted as a function of the temperature T. As can 
be seen, the Boltzmann approach gives correct orders of 
magnitude in the temperature range and for the CNT 
index considered here. It however overestimates the con- 
ductivity, and the error increases with the length of the 
system and the temperature. This can be explained 
by the fact that the Boltzmann approach is not precise 
enough to capture the real decay of the transmission func- 
tion. In particular, the decrease of the transmission is not 
fast enough, which is consistent with the results of Fig. [S] 



However, as already pointed out in Section IV Bl it is 
expected that the Boltzmann treatment becomes more 
accurate as the number of interacting modes increases. 
It is therefore likely that the agreement between the 
conductance computed using Green's function techniques 
and conductances computed from the Boltzmann trans- 
missions will be better for tubes with larger indices or 
multi-walled tubes (or for models treating the out-of- 
plane modes more precisely). This is indeed true, as can 
be seen from the results presented in Fig. [T31 where the 
error on the conductivity is smaller for CNTs of larger 
diameters (To give some orders of magnitude, the diam- 
eters of (5,5) and (10,10) CNTs are respectively 0.69 nm 
and 1.38 nm). 
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FIG. 13: Error on the the conductivity computed with the 
profile obtained from the Boltzmann approach. The reference 
value is computed from the profiles obtained with the Green's 
function approach (decreasing lengths from top to down). 
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VI. CONCLUSION 

We studied the thermal transport of isotopicaUy disor- 
dered harmonic CNTs using Green's function and Boltz- 
mann treatments. The Green's function method is im- 
plemented using an algorithm, whose computational cost 
scales linearly with the system size. It is not restricted 
to isotopic disorder, and can be applied to other kind 
of defects which scatter phonons elasticallyi^ In view of 
the theoretical results for ID systems, we also expect a 
divergence of the thermal conductivity with the system 
size in those cases. 

Our numerical results show that 

(i) these systems share some common features with 
the thermal transport in isotopicaUy disordered har- 
monic one dimensional atom chains, in particular a 
power law divergence of the thermal conductivity 
K ~ L". Therefore, Fourier's law is not valid. For 
tubes of experimental length, the exponent of the 
power law divergence a ~ 0.6 at room temperature. 
For longer tubes, the exponent may well be differ- 
ent since a = 1/2 in the asymptotic regime where 
only the acoustic branches have a non-zero trans- 
mission. Experimental results on the length depen- 
dence of disordered CNT conductivities, measuring 
divergences a ~ 0.4 — 0.6 for boron nitride tubes 
(with a large fraction of mass disorder), have been 
published in Ref . @ 

(ii) the thermal conductivity decreases monotonically 
with the temperature, and is independent of the 
CNT diameter. We showed that there is a dramatic 
reduction of the thermal conductance for systems of 
experimental sizes (roughly 80% at room tempera- 
ture), when a large fraction of isotopic disorder is 
introduced. This is in accordance with experimental 
measurements of the effect of isotope disordered in 
boron nitride nanotubes, which demonstrated dra- 
matic changes in the conductances (a conductance 
enhancement by 50% for purified materials) 

(iii) a Boltzmann description of the thermal transport 
gives correct results for the thermal conductance 
even in the presence of Anderson localization. This 
is particularly interesting since the computation 
of the transmission using Boltzmann's equation is 



much less computationally expensive, so that larger 
systems (such as multi-walled CNTs or boron ni- 
tride nanotubes, or single- walled CNTs with a large 
diameter) may be studied with this method. This 
shows also that Anderson localization, that is not 
accounted for in the Boltzmann approach, does 
not have a clear signature in thermal transport for 
CNTs of experimental sizes; 

(iv) the results discussed in the previous points are ro- 
bust with respect to the model used to describe the 
phonon dispersion in the nanotubc. In particular, 
the behaviors predicted by a simple empirical model 
which neglects the tube curvature are very similar to 
those obtained with ab-initio force constants, com- 
puted with density functional theory, and consider- 
ing the tube curvature. 
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APPENDIX A: COMPUTATION OF THE 
TRANSMISSION WITHIN THE GREEN'S 
FUNCTION FORMALISM 

In order to simplify the notations, the number of layers 
in the disordered region A^iaycrs will be denoted by N in 
this section. Besides, the parameter rj is always kept pos- 
itive in numerical computations. We indicate explicitly 
the dependence on this parameter in this section. 

The {i,j) block component of some matrix Ai is de- 
noted by It is easily seen that the only non zero 
block components of the self-energies are respectively 
[S^((jj)]i.i and [S^(tj)]Arjv. These coefficients are com- 
puted using standard decimation techniques.— >^i^i^ 

The transmission T{uj) can be evaluated using only the 



block [G^ (tj)]i^Ar of the full matrix Glj (to) 



J 



T{cj) = 4Tr [lm([S^(^)]i,i) [G^,,(c.)]i,^ Im([S^(c.)]^,^) {[Gl.icj)],^^^ 



where the trace is taken over a space of dimension 3iVat x ing 

3iVat. 

The submatrix [G^^^g{uj)]i^N can be computed by solv- 












/ 
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with 



We use a method to compute hnear scahngly the last 
column of the inverse of a block tridiagonal matrix, 
based on gaussian elimination. Each [G''ly^{uj)]i^N is 
expressed in terms of [G^ (a;)]i+i,jv as [G^ {uj)]i^M ~ 



F,[G2^^{oj)Ui,N- Denoting by B 



IT] - As 



T.l{u) - and /, = [G,^y,(w)],,w, it holds, for 

2<i<N -I, 

Bis-ifi^i + Bjjfi + Bi^i^ifi+i =0, fi = Fifi+i. 

It is then easily shown that 



F, 



while /tv = {Bn,n + -Bjv,jv-i^jv-i) ^, so that 
[G7ys{^)]i,N /i = F1F2 . . .Fn-iJn is recovered in 
0{Ncfi) operations. The corresponding algorithm there- 
fore has a cost O(Mayers-^at) (^^^ Algorithm [T]) . 



Iterative computation of the {1,N) block of the 
effective resolvent 

Algorithm 1. Fix lo^ > and rj > 0, and denote B — 
uj" + iri ~ Asys — T,^{lu) — E^(tj). Set 

Fi = -B^}Bi,2, 9 = Fi. 

Then, 

(1) For i = 2, . . . , N — 1, compute 

Fi = ~{Bi^i + Bi^i-iFi-i) ^ Bi^i+i, 

and replace g by gFi; 

(2) Set [G'',y,{uj)]i,N = g{BN,N + Bn,n-iFn-i)-\ 



An alternative linear scaling algorithm to compute the 
transmission is proposed in Ref. 1681 . This algorithm also 
allows to compute the density of states in O(A^iayorsA'^fj) 
operations, but gives wrong transmissions for non small 
values of rj (in the results presented in this work, this 
occurs for even for values ry/w^ax = 10~^ for one dimen- 
sional chains long enough or when the variance of the 
mass disorder is large enough). 



APPENDIX B: NUMERICAL 
IMPLEMENTATION OF THE BOLTZMANN 
EQUATION 

To simplify the notations, we drop the variable lo in 
the phonon distributions (since it merely indexes inde- 
pendent Boltzmann equations). To solve (|14p . it is con- 
venient first to symmetrize the problem by introducing 



fj,a = VjUj^a, so that (jl4p can be recast as 

(Bl) 

where the matrix A is symmetric. More precisely, 



W, 



(j,'T),(j',cr') 



when (/, a') 7^ (j, a) and 



A 



Wi 



The boundary conditions for (jBip are obtained from the 
boundary conditions for ()14p : 



/,,+ (0,i)-i;„ /,,-(L,i) = 0. 

We denote / = (/i,+, . . . , /Ar,+, . . . , /jv,-)- Recall 
that 2iV is the number of conducting channels at this 
pulsation. The stationary solution of (|B1[) is such that 
!{x) = exp(a;^)/(0), and thus /(O) = exp(-L^) /(L), 
where we dropped the argument t. This expression has 
a meaning since A is real symmetric hence defines a self- 
adjoint operator. Expressing exp(— LA) as a 2 x 2 matrix 
where the entries are N x N matrices: 



exp(-LA)- f B22(S 



it holds Bn{L) ■ (/i,+ (L), . . . , fN.+{L)Y ^{vi,..., vnY- 
Defining the N x N matrix B{L) ~ B^^{L) when this 
inverse exists, the transmission is finally 



T ^ N 



i=i 



N 

E^ 



TV 



1=1 j=i 



N 

E^ 



In practice, it is not possible to use this formula, since 
the matrix Bii(L) becomes singular when L is large. 
In those cases, we compute a stationary solution of the 
Boltzmann equation (|B1[) by integrating in time the cor- 
responding partial differential equation (PDE). The dis- 
cretization is done on a reference geometry (thanks to 
the hyperbolic scaling of the transport part, it is possible 
to map a system of size L to a system of size unity, upon 
rescaling time and the scattering terms by a factor L). 
The numerical scheme uses a splitting of the dynamics 
into a transport part (which is implemented with a sim- 
ple upwind scheme^ since the sign of the velocities are 
known and are constant), and a coUisional part (which 
is integrated exactly using a matrix exponential). The 
stability conditions on the time step (the so-called CFL 
conditions in the literature on systems of conservation 
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laws^) arise therefore only from the transport part of 
the equation, and are easy to state since the geometry is 
kept fixed. Wc have checked that the algebraic method 
using the inversion of matrix exponentials and the com- 
putation of the steady state through the time integration 
of the PDE give the same results when the condition 
number of the matrix Bn is not too large. Another con- 
sistency check is the fact that the computed transmission 
is almost constant in space (up to the convergence error). 

In the simulation results presented in this work, we 
have used = 200 discretization points on the refer- 
ence interval [0,1], and a time step At = 1.25 x 10"'^ 
for the time integration of (|B1[) . We have checked the 
convergence of the results with respect to Nx and At. 

The scattering matrices (fTS]) have been obtained by 
considering a discrete phonon spectrum computed with 
10** fc-points in the first Brillouin zone, and by using win- 
dows of width 3.2 cm""'^ centered around the pulsations 
Lu of interest to list all the branches at this pulsation. An 
approximation of the phonon velocity is obtained through 
finite differences. Some scattering matrices are discarded 
in the end when they contain too large values (which 
happens when the velocity of a branch vanishes). This 
represents less than 1% of the cases for the profiles com- 
puted. Such a procedure cannot however be used when 
the CNT index increases since more and more branches 
start with a slope zero, and refined strategies should then 
be employed for larger systems. 



n_|_(x) — n^{x) 



where the dependence on uj is not written out explicitly, 
and 



^-1 ^ ^(^^-1 ^ W+A^) _ ^Var(m) uj\k) 



From (|C1[) . dx{n+ — ?i_) =0. Using the boundary condi- 
tions n.-|_(0) = 1, n-[L) ~ 0, and denoting by n^{L) = T 
and n_ (0) = R the proportions of transmitted and re- 
flected waves respectively, it follows (ri-^. — n^){x) = 
(n_|. — n_)(0) ~ (n-|_ — n^){L) = T. Plugging this re- 
suh in (ICTI). 



dxn+{x) 



T 



so that 



The transmission is then computed at x ^ L using the 
above expression for 7i+: 



APPENDIX C: DERIVATION OF THE 
BOLTZMANN FORMULA ([T9|) 

For one dimensional chains, there is a single conducting 
branch parametrized by k in the first Brillouin zone: 



Lu{kf = 2(1 - cos(fca)). 



<k < -, 

a a 



(CI) 



where a denotes the lattice spacing in this appendix (In 
section ITVl a = 1). The steady-state Boltzmann equation 
at some frequency to reduces to 



dxn+{x) 



n_(x) — n^{x) 
I 



Tl{u;) 



l(u}) + L' 



To obtain (|19p . it remains to precise the value of I ^ 
l{Ld)~^ . From (|C1|) . the phonon velocity is 



uj'{k) 



as\Q.{ka) 



so that the result follows using sin(fca)^ ~ 1 — cos(fca)^ 

1- (l-tjV2)2 =Cj2 (1-^^2/4). 
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